Stable reference gene selection for Ophiocordyceps sinensis gene expression studies under different developmental stages and light-induced conditions

The molecular mechanism of Chinese cordyceps formation has received a substantial amount of attention because of its usage as traditional Chinese medicine. The formation process of Chinese cordyceps includes two parts: asexual proliferation (Ophiocordyceps sinensis proliferates in the hemolymph of Thitarodes armoricanus larvae) and sexual development (formation and development of fruiting bodies). Therefore, validation of reference genes under different development stages and experimental conditions is crucial for RT-qPCR analysis. However, there is no report on stable reference genes at the development stage of O. sinensis fruiting body. In this study, 10 candidate reference genes, Actin, Cox5, Tef1, Ubi, 18s, Gpd, Rpb1, Try, Tub1 and Tub2, were selected and calculated their expression stability using four methods: geNorm, NormFinder, BestKeeper, and Comparative △Ct. After comprehensive analysis of the results of these four methods with RefFinder, we determined that the most stable reference genes during asexual reproduction of O. sinensis were Tef1 and Tub1, while the most stable reference genes during fruiting body development were Tyr and Cox5, and the most stable reference genes under light-induced conditions were Tyr and Tef1. Our study provides a guidance for reference genes selections at different proliferation processes with light stress of O. sinensis, and represents a foundation for studying the molecular mechanism of Chinese cordyceps formation.


Introduction
Chinese cordyceps, also named "Dong Chong Xia Cao" is a traditional Chinese medicine. It is an insect-fungus complex formed by an entomopathogenic fungus, Ophiocordyceps sinensis, that parasitizes the larvae of the insect Thitarodes armoricanus. Increasing in-depth research has revealed that Chinese cordyceps shows many pharmacological anti-oxidative [1], antiinflammatory [2], hypoglycemic [3], and immunomodulatory [4] actions. Chinese cordyceps is endemic to the Qinghai-Tibet Plateau and mainly occurs at altitudes between 3000-5000 m.

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0284486 April 20, 2023 1 / 13 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Extreme growth environments limit wild Chinese cordyceps yields. Therefore, researchers have begun to study artificial cultivation of Chinese cordyceps in low altitude areas and have made many breakthroughs [5,6]. In the low-altitude laboratory, the researchers were able to complete the inoculation of O. sinensis in the larval stage of the T. armoricanus, and induce larval death and mummification to grow fruiting bodies [7]. However, low mummification rates and inconsistent fruiting body development are still major problems that affect artificial cultivation [8]. Consequently, many researchers have focused on fruiting body formation [7,9] and the influence of environmental factors (such as temperature and light) on fruiting body development [10,11]. Through transcriptome and proteome analyses, many differentially expressed genes were identified in different stages of fruiting body development [12]. These genes are involved in MAPK signal transduction, glucose metabolism [13], amino acid metabolism [14], lipid metabolism, and other pathways [15,16]. Simultaneously, it was also found that O. sinensis is sensitive to light [11], and the mycelia synthesize melanin under light conditions (especially under UV irradiation) [17]. A large amount of omics data provided many new targets for O. sinensis artificial cultivation. However, no stable reference gene has been reported, which is very important for studying the molecular mechanism. Studying the key genes and molecular mechanisms in the process of O. sinensis infection, asexual proliferation, sexual development, and synthesis of active substances will become the research focus of artificial cultivation of Chinese cordyceps. During the formation of Chinese cordyceps, we usually focus on the formation of fruiting bodies and the effect of light on O. sinensis.
Reverse transcription quantitative PCR (RT-qPCR) is a widely used and authoritative method to evaluate gene expression. However, some variable factors may affect the accurate quantification of gene expression [18,19], including the differences in materials, RNA extraction process, reverse transcription efficiency, and primer design. Therefore, stably expressed reference gene is needed to correct experimental errors [20]. It is generally believed that housekeeping genes are stably expressed and used for standardized correction of target genes. However, many studies have proved that the expression of common housekeeping genes will change due to different species and experimental treatments, and there were no konwn stable reference gene applicable to all conditions. Thus selections of a collection of reference genes is necessary for particular experimental models. The formation process of Chinese cordyceps involves the asexual and sexual multiplication stage of O. sinensis. In the stage of asexual proliferation, O. sinensis proliferates in the host hemolymph, and there will be morphological transformation between blastospores and hyphae [12]. In the stage of sexual development (fruiting bodies formation and maturity), similar to other edible fungi, light is usually used to regulate the phenotype of fruiting bodies. In this process, there is no report on the stable reference gene of O. sinensis.
In order to fill this gap, we explored the conditions of liquid fermentation and fruit-body induction, and observed and prepared O. sinensis samples at all stages of development as comprehensively as possible. Liquid culture were used to simulate the process of asexual proliferation, and through this method, blastospores and hyphae (or mycelial pellets) of O. sinensis under different light conditions were collected. In the process of sexual reproduction, we prepared a batch of T. armoricanus larvae infected with O. sinensis. Under the condition of artificial induction, the hyphal bodies in hemolymph and the fruiting bodies at different stages of development (form formation to maturity) were collected. By using the methods of geNorm [21], NormFinder [22], BestKeeper [23], and Comparative 4C t methods [24], the expression stabilities of ten common reference genes was analyzed.We screened the stable reference genes of O. sinensis in the asexual proliferation process, the fruiting body development process and the light-induced process. Our results will provide guidelines for reference gene selection of RT-qPCR in O. sinensis gene expression analysis

Strains cultivation and sample collection
For the asexual samples, O. sinensis was activated and cultivated with OS1 liquid medium (20% potato infusion, 2% glucose, 1% yeast extract, 0.5% peptone, 0.1% KH 2 PO 4 , 0.025% MgSO 4 ), and shaken at 16˚C and 100 rpm for 40 d. The blastospore suspension was collected as seed solution with lens wiping paper, and inoculated in OS1 and 1/4 SDY liquid medium (1% glucose, 0.5% yeast extract, 0.25% peptone), and then shaken at 16˚C and 120 rpm. Light and dark treatments were performed on the samples cultivated with OS1 liquid medium. The blastospore and mycelium samples were harvested in 1/4 SDY liquid medium at 30 d and 45 d, respectively. The loose and dense mycelium balls were harvested in OS1 liquid medium at 30 d and 70 d, respectively.
The larvae infected with O. sinensis and the fruiting bodies were reared in a low-altitude laboratory (Jiangjin, Chongqing). The hyphal body samples were collected from the hemolymph of larvae inoculated with O. sinensis. Mummified larvae were buried in culture medium for moisture retention during the fruiting body formation and development. Light (light: dark = 12 h: 12 h) and dark (24h dark) treatments were performed on the samples after fruiting body formation.

RNA isolation and cDNA synthesis
The samples were ground in liquid nitrogen, and 200 mg samples were used to isolate total RNA with the Total RNA Kit II (R6934-02, Omega, Norcross, Georgia USA). Then, 2 μg total RNA of each sample was used for reverse transcription synthesis of cDNA with Hifair1 Ⅲ 1st Strand cDNA Synthesis SuperMix for qPCR (YEASEN, Shanghai, China) and stored at −20˚C after 10-fold dilution. qPCR primers were designed by the NCBI Primer-BLAST tool (https:// www.ncbi.nlm.nih.gov/tools/primer-blast/).

Candidate reference gene selection and primer design
According to similar studies in other fungi [18,19,25,26], ten housekeeping genes were selected as candidate reference gene targets, including actin protein (Actin), cytochrome c oxidase polypeptide V (Cox5), translation elongation factor EF-1 alpha (Tef1), polyubiquitin binding protein (Ubi), 18s ribosomal (18s), glyceraldehyde-3-phosphate dehydrogenase (Gpd), polymerase II large subunit (Rpb1), tryptophan synthase (Try), and tubulin beta chain (Tub1, Tub2); the tubulin beta chain has two copies in O. sinensis that we named Tub1 and Tub2. The candidate genes covered multiple processes such as transcription and translation (Rpb1 and Tef1), protein synthesis and degradation (18s and Ubi), cytoskeletal structure (Actin, Tub1, and Tub2), and metabolism (Cox5, Gpd, and Try). Details of reference genes and primers are shown in Table 1. Specificity of primer pairs for each candidate gene was confirmed using melting curve analysis (S1 Fig).

RT-qPCR
The samples were divided into three groups: asexual proliferation (all liquid fermentation samples), sexual proliferation (samples during fruiting body development), and light treatment (light-and dark-treated samples during liquid fermentation and fruiting body development). Quantitative PCR was used to detect the expression of ten reference genes in each sample with LightCycler1 96 (Roche, Indianapolis, IN, USA). The PCR mix (10 μL) contained 5 μL Hieff 1 qPCR SYBR Green Master Mix (YEASEN, Shanghai, China), 0.5 μL per primer, 1 μL 10-fold diluted cDNA, and 3 μL ddH 2 O. The PCR conditions were as follows: 95˚C for 5 min, followed by 40 cycles of 95˚C for 10 s and 60˚C for 30 s. After the final amplification cycle, the melting curve analysis was performed as follows: 65˚C for 60 s, 65-95˚C with 0.5˚C increments, and 97˚C for 10 s.

Statistical analysis
The expression stability of ten candidate reference genes was analyzed using geNorm, Norm-Finder, BestKeeper, and Comparative 4C t methods. GeNorm program converts cycle number (Ct) into relative gene expression (2 ΔCt ), evaluated by calculating the stable value (M) of candidate reference gene expression, and finally, two or more reference gene combinations were provided by the analysis results of pairwise variation (Vn/Vn+1) [21]. NormFinder is also sorts the stability of reference genes according to M, and only one stable reference gene can be provided [22]. BestKeeper evaluated the stability of reference gene expression according to standard deviation (SD) and coefficient of variation (CV), the more stable genes have lower SD values [23].4C t method is to compare the 4Ct values of paired genes in different samples. If the 4C t values of these two genes remain unchanged in all samples, it is considered that these two genes are stably expressed in these samples; Otherwise, the expression of one or two genes is considered unstable. The final gene stability is sorted according to the change of 4C t value [24]. The overall final ranking was calculated by RefFinder(http://blooge.cn/RefFinder/) [27], which assigns the appropriate weights to individual genes based on the ranking of each program and calculates the geometric mean of their weights for the overall final ranking. The data analysis was performed using LightCycler1 96 Software 1.1 (Roche), and visualization using Origin. Each of the above experiments was independently repeated at least 3 times.

Sample collection
Through liquid fermentation, we collected O. sinensis samples with five morphologies (Fig 1A). They are blastospores, hyphae, loose hyphal spheres, and compact hyphal spheres under light  Fig 1B-f, light: Fig  1B-g), 18 d (dark: Fig 1B-h, light: Fig 1B-i), 32 d (dark: Fig 1B-j, light: Fig 1B-k), 47 d (dark: Fig 1B-l, light: Fig 1B-m), and 62 d (dark: Fig 1B -n, light: Fig 1B- and dark conditions. In the early stage of fermentation, 5-10 μm fusiform blastospores were proliferated by budding (Fig 1A-a red arrows). After continued cultivation for 45 d, most of the blastospores germinated into hyphae (Fig 1A-b). In OS1 medium, the hyphae aggregated into mycelial pellets (Fig 1A-c). The mycelial pellets become dense and smooth after 70 days of culture, and appear white (Fig 1A-d) and black (Fig 1A-e) respectively under the induction of dark and light conditions. The hyphal body samples were harvested from larval hemolymph (Fig 1B-a) and were observed as fusiform-like blastospores. The hyphal body length (10-15 μm) was longer than that of the blastospores (5-10 μm). The mummified larvae were embedded in the culture medium to induce fruiting body formation. Initially, O. sinensis only proliferates in the host, and the mycelium does not penetrate the host epidermis so that the head shell of the mummified larvae (Fig 1B-b) did not change significantly. However, the fungus preferentially penetrated the larval epidermis at the head shell with proliferation and form a 'Y-type' hyphal structure (Fig 1B-c). Subsequently, the 'Y-type' hyphal structure swelled and bulged (Fig 1Bd), and finally developed into the primordium (Fig 1B-e). The results of the mummified larvae when primordia exposed to light and dark treatments showed that there was no significant difference in fruiting body phenotypes between dark (Fig 1B-f) and light (Fig 1B-g) treatments in the early stage of fruiting body development. When the fruiting body developed to approximately 2 cm, the color of the top of the fruiting body became yellow-brown with light treatment (Fig 1B-i) but was still white under dark treatment. Upon continued cultivation of the fruiting body to 7 cm, the color of fruiting body was yellow-brown under light treatment ( Fig  1B-k), whereas the color only got darker at the top and base under dark treatment (Fig 1B-j). After 40 days of fruiting body development, the middle and upper part of the light-treated fruiting bodies formed bulges (Fig 1B-m) and finally developed into asci (Fig 1B-o), whereas the dark-treated fruiting bodies still had a smooth surface (Fig 1B-i, n).

Expression profiling of candidate genes
The expression levels of 10 genes were calculated by quantification cycles and shown with box plots (Fig 2). The samples were divided into three groups (asexual proliferation, sexual proliferation, and light treatment), and the candidate genes had the same expression trend in different groups. Except for the Cq value of 18s distributed around 10, the expression levels of other genes were in the range of 20-30 indicating that the expression level of 18s was significantly higher than that of other genes. In the sexual proliferation and light treatment groups,

geNorm analysis
The average expression stability (M value) of the ten candidates was calculated with the GEN-ORM software. The M value is defined as the average pairwise variation of a particular gene with all other reference genes within a given group of cDNA samples. where a low value represents stable and a high value represents unstable expression. A low stability value (M) indicates high expression stability and the geNorm software defaults that M = 1.5 is the threshold of gene stability. According to the calculation results, 7 genes in asexual proliferation samples were within the acceptable range (Table 2), and the most stable reference genes are Tub1 / Tyr (M = 0.153). Similarly, there are 6 genes within the acceptable range in the fruit body development samples (Table 3) and the light stress samples (Table 4), and the most stable reference genes are Cox5/ Gpd (M = 0.619) and Tyr / Tef1 (M = 0.579) respectively.
Further, the minimal number of genes for qRT-PCR normalization were determine by estimated pairwise variation (Vn/Vn+1) in geNorm. It is generally considered that Vn/Vn +1 < 0.15 is the cut-off threshold, and no additional reference gene is required. Therefore, no more than 2 reference genes are required in ll three groups of samples (0.007 in asexual proliferation samples, 0.012 in fruit body development samples and 0.011 in light stress samples) to normalization the qRT-PCR results (Fig 3).

NormFinder analysis
NormFinder calculate the stability based on an ANOVA-based algorithm, and also considers the best two genes with minimal combined intra-and intergroup expression variation. In 1 The stability value calculated by geNorm is the standard deviation of the expression variation of the two genes. A low stability value indicates high expression stability. 2 Comparative ΔC t method ranks candidate gene stability according to the repeatability of gene expression differences. A low Average of STDEV indicates high expression stability. 3 The stability value calculated by NormFinder. A high stability value indicates low expression stability. 4 BestKeeper analyzes expression stability by calculating cycle threshold (Ct) data variation. A low SD (standard deviation) and CV (Coefficient of Variance) value indicates high expression stability. 5 Comprehensive stability value assigns the appropriate weights to individual genes based on the ranking of each program and calculates the geometric mean of their weights for the overall final ranking. 6 Ranking assigned to each gene according to its comprehensive stability value.  (Table 3). In light stress samples, the order of stability (most stable to least stable) is (Table 4). In the calculation results of NormFinder and geNorm, the most and least stable reference genes are basically the same.

BestKeeper analysis
BestKeeper evaluates the stability of reference genes by calculating the standard deviation (SD) and the coefficient of variation (CV) of Cq values. It is generally believed that the reference genes with SD < 1 are stable, and the smaller the CV indicates the higher stability. As the result, 7 genes (Tef1, Actin, Tub2, Tyr, Tub1, Rpb, 18s) in asexual proliferation samples (Table 2), 2 genes (Rpb and 18s) in fruiting body development samples (Table 3) and 3 genes (Rpb, Tub1, 18s) in light stress samples (Table 4) are stable. The results of BestKeeper are different from NormFinder and geNorm, so it is necessary to evaluate the results of multiple software. The footer is the same as Table 2 https://doi.org/10.1371/journal.pone.0284486.t003 The footer is the same as Table 2 https://doi.org/10.1371/journal.pone.0284486.t004

Comparative ΔCt analysis
The comparative ΔCt method evaluate reference gene stability by the repeatability of gene expression differences. In asexual proliferation samples (  (Table 3) and light stress samples samples (Table 4).

Comprehensive ranking order
The reference genes stability calculated by different methods is slightly different, so RefFinder was used to comprehensively calculate the overall final ranking. RefFinder assigns the appropriate weights to individual genes based on the ranking of each program and calculates the geometric mean of their weights for the overall final ranking. we recommended that the most stable reference genes during asexual reproduction of O. sinensis are Tef1 and Tub1 (Table 2), the most stable reference genes during fruiting body development were Tyr and Cox5 (Table 3), and the most stable reference genes under light-induced conditions were Tyr and Tef1 (Table 4).

Discussion
RT-qPCR is an important method for studying gene expression that can help eliminate experimental errors, and selecting a stably expressed reference gene for RT-qPCR analysis is necessary [26]. In this study, the asexual proliferation process, fruiting body development process, and effects of light on O. sinensis were recorded in detail. Additionally, we assessed the stability of 10 candidate reference genes under three different experimental conditions. These results provide a reference for studying the gene expression of O. sinensis under different conditions. Through this study, we found that the reference genes stably expressed by O. sinensis are not the same under different proliferation forms and treatment conditions, which further demonstrated the importance of selecting appropriate reference genes. The formation process of Chinese cordyceps can be divided into two stages. One is the proliferation stage of O. sinensis in the hemolymph of T. armoricanus larvae. It is generally believed that, in this stage, O. sinensis mainly propagates asexually in the form of hyphal bodies [28]. At the subsequent stage, the hyphal bodies are induced by certain signal substances to transform by hyphae proliferation [29]. Our results revealed that, as the content of O. sinensis mycelium in the host hemolymph increased, the mycelium began clustering and caused the host larvae died and mummified. The second stage is the formation of fruiting bodies [12]. It is generally believed that, similar to mushrooms [30], the formation of fruiting bodies of O. sinensis is also the result of sexual mating, and asci will be formed in the middle parts of the fruiting bodies at the later stage of development [12]. Although the molecular mechanism of Chinese cordyceps formation has not been completely analyzed, the dimorphic transition (from blastospores or hyphal bodies to mycelial pellets) [31] and the developmental regulation of fruiting bodies have attracted the attention of many scholars.
In this study, the asexual proliferation process of O. sinensis was simulated by liquid fermentation. The formation of blastospores (similar to hyphal bodies), mycelia, and mycelial pellets was more accurately reproduced. Using different types of samples, we screened the stably expressed reference genes Tef1 and Tub1 during asexual proliferation. Although it could not be proven that the proliferation process of O. sinensis in liquid fermentation and host hemolymph was the same, the similar morphological changes can also provide a reference for studying the dimorphic transition of O. sinensis [29]. Our findings provide a tool for future research on morphogenesis and morphological changes of O. sinensis.
The formation and development of fungal fruiting bodies involve complex gene regulatory networks, and there are great differences among different species. Therefore, different experimental models need to re-evaluate the stability of internal reference genes. In Cordyceps militaris, the best reference genes for inducing fruiting bodies in wheat medium and pupae are gpd/rpb1 and rpb1/try, respectively [32]. However, during the fruiting bodies development of Flammulina velutipes, Botrytis cinerea and Agaricus bisporus, the stable reference genes are gapdh, tubulin and EF1-α, respectively [33][34][35]. In the research related to the development of O. sinensis fruiting body, no research has evaluated the best reference genes yet. In this study, we collected samples from the mummified larva until the fruiting body matured and formed asci. By using fruiting body samples from various developmental stages, we determined that Tyr and Cox5 represented stably expressed reference genes during the fruiting body development process.
During the process of asexual proliferation and fruiting body development, we found that light can affect O. sinensis color and ascus formation. Fungi can respond to light in various ways. For examples, fungi can respond to different wavelengths of light through photoreceptor proteins, and can also regulate the expression of downstream genes through light-induced intracellular ROS and NOS [36,37]. In C. militaris, researchers usually use Rbp1 as an internal reference gene when studying the regulation of light in fungal gene expression [32]. However, in this study, we found that Rbp1 is not a stable reference gene for O. sinensis response to light. In a previous study, the researchers used 18s as a reference gene to explore the expression pattern of the blue light receptor OsWC1 in O. Sinensis [38]. It is generally believed that 18s expression is stable and not easily affected by other conditions, but the expression level of 18s is very high compared with other genes, and it is likely to exceed the confidence interval of the standard curve [39,40]. In this study, we determined that the reference genes Tyr and Tef1 are more stable than 18s and Rbp1, which can be used for follow-up studies of O. sinensis responding to light.

Conclusions
This study comprehensively described and collected samples of various phenotypes of sexual and asexual proliferation of O. sinensis, including blast spores, hyphae and hyphal spheres in liquid fermentation, hyphal body in the host hemolymph and fruiting bodies at various development stages. At the same time, the influence of light on each stage was studied. Moreover, three groups of samples (asexual proliferation, fruiting body formation and development, and response to light stress) have been screened for stable reference genes that are suitable for different experimental models through four methods (geNorm, Comparative 4C t , NormFinder, and BestKeeper). This is the first time to screen reference genes comprehensively in O. sinensis.
This study provides a reliable tool for further studying the gene expressions in the process of O. sinensis morphological transformation, fruiting body development, and response to light stress.